Varroa EAG response to various essential oils
Nurit Eliash
3/16/2022
Load libraries
library(dplyr)
library(tidyverse)
library(ggplot2)
library(ggpubr)
library(ggthemes)
library(stats)
library(moments)
library(grid)
library(formattable)
library(gridExtra)
library(ggsignif)
library(hrbrthemes)
library(plotrix)
library(rstatix)
library(car)
library(plotly)
library(postHoc)Load data
data <- read.csv("/Users/nuriteliash/Documents/GitHub/Sofia-varroa-project/EAG_data_all.csv")EAG stat analysis
the tested effect (x) is the Essential oil “EO” dissolved in acetone, the “response” (y) is the electrophysiological response of the leg measured in mV. Each leg was stimulated by different stimuli in the following order:Air > Air > Air > Acetone (solvent) > 0.1 > 0.25 > 0.5 > 1 > 2.5 > 5 > Acetone (solvent) > Air
the response was normalized…. [please complete in here, by copy-pasting from the method section:)]
to test the difference in response amplitude to the different stimuli, we will use One way ANOVA. then a post-hoc Tukey test, to see which of the stimuli are significantly differnet. we followed this tutorial for ANOVA in R
Data summary
# our data contains:
str(data)## 'data.frame': 664 obs. of 4 variables:
## $ leg : Factor w/ 83 levels "VF1","VF10","VF11",..: 1 1 1 1 1 1 1 1 12 12 ...
## $ stimuli : Factor w/ 8 levels "0.1","0.25","0.5",..: 8 1 2 3 4 5 6 7 8 1 ...
## $ EO : Factor w/ 9 levels "EO4309","EO4444",..: 2 2 2 2 2 2 2 2 2 2 ...
## $ response: num 109.5 118.5 86.6 131.9 108.5 ...table(data$stimuli, data$EO)##
## EO4309 EO4444 EO4891 EO4899 EO5068 EO5663 EO5714 EO5749 EO5763
## 0.1 8 10 10 10 11 9 9 9 7
## 0.25 8 10 10 10 11 9 9 9 7
## 0.5 8 10 10 10 11 9 9 9 7
## 1 8 10 10 10 11 9 9 9 7
## 2.5 8 10 10 10 11 9 9 9 7
## 5 8 10 10 10 11 9 9 9 7
## acet_after 8 10 10 10 11 9 9 9 7
## acet_before 8 10 10 10 11 9 9 9 7We have tested a total of 9 essential oils, using 83 legs. Each essential oil was tested on at least 7 different legs.
Prior to analysis, we gonna check the ANOVA assumptions:
(1) Outliers
# plot it to detect outliers by specific leg
# first sort the order of stimuli:
data$stimuli <- factor(data$stimuli,levels = c("acet_before", "0.1", "0.25", "0.5", "1", "2.5", "5","acet_after"))
box <- ggplot(data, aes(x = stimuli, y = response)) +
geom_boxplot(aes(colour = EO)) +
facet_wrap( ~ EO) +
theme_linedraw() +
ggtitle("Varroa foreleg response to different essential oils") +
xlab("Stimuli amount (microgram)") +
ylab("Normalized response (%)") +
theme(axis.text.x=element_text(angle=45, hjust=1))
ggplotly(box, tooltip = "leg")Remove outliered legs: “VF71”, “VF1”, “VF2”, “VF4”, “VF77”, “VF82”, “VF24”, “VF26”, “VF37”, “VF39”, “VF46”, “VF49”, “VF17”, “VF19”, “VF50”, “VF62”
data <- data %>%
dplyr::filter(!leg %in% c("VF71", "VF1", "VF2", "VF4", "VF77", "VF82", "VF24", "VF26", "VF37", "VF39", "VF46", "VF49", "VF17", "VF19", "VF50", "VF62"))After excluding the outliers, we proceed with the rest of the tests:
(2) Normality:
#the dependent variable should be approximately normally distributed in each cell of the design. This can be checked using the Shapiro-Wilk normality test (shapiro_test()
normality <- data %>%
group_by(stimuli,EO) %>%
shapiro_test(response)
normality %>%
dplyr::filter(p<0.05)## # A tibble: 4 x 5
## stimuli EO variable statistic p
## <fct> <fct> <chr> <dbl> <dbl>
## 1 acet_before EO4309 response 0.761 0.0163
## 2 0.1 EO5663 response 0.792 0.0344
## 3 1 EO5068 response 0.783 0.0130
## 4 acet_after EO4891 response 0.738 0.00607For the above four EO x stimuli combinations the response dont distribute normally (p<0.05)
(3) Homogneity of variance
# Build the linear model
model <- lapply(split(data, data$EO), function(i){
lm(response ~ stimuli, data = i)
})
# now you can Create a QQ plot of residuals of each EO, for example, essential oil EO4309:
ggqqplot(residuals(model$EO4309))
plot(model$EO4309, 1)
After checking the 3 assumptions, we decided to go for a non-parametric test, to test the significant difference of each stimuli concentration from the solvent.
we used a non-parametric post hoc test for multiple comparisons.
Post hoc non-parametric test for all pairwise comparisons, with Benjamini-Hochberg p-value correction
Post hoc Kruskal-Wallis rank sum test
KW <- lapply(split(data, data$EO), function(i){
posthocKW(i$response, i$stimuli, padjust = "BH")
})
KW$EO4309$PvaluesMatrix## acet_before 0.1 0.25 0.5 1 2.5
## acet_before NA NA NA NA NA NA
## 0.1 0.9490597 NA NA NA NA NA
## 0.25 0.9490597 0.9490597 NA NA NA NA
## 0.5 0.9490597 0.9490597 0.9490597 NA NA NA
## 1 0.9490597 0.9490597 0.9490597 0.9490597 NA NA
## 2.5 0.9490597 0.9490597 0.9490597 0.9490597 0.9490597 NA
## 5 0.9490597 0.9490597 0.9490597 0.9490597 0.9490597 0.9490597
## acet_after 0.9490597 0.9490597 0.9490597 0.9490597 0.9490597 0.9490597
## 5 acet_after
## acet_before NA NA
## 0.1 NA NA
## 0.25 NA NA
## 0.5 NA NA
## 1 NA NA
## 2.5 NA NA
## 5 NA NA
## acet_after 0.9490597 NAKW$EO4444$PvaluesMatrix## acet_before 0.1 0.25 0.5 1 2.5
## acet_before NA NA NA NA NA NA
## 0.1 0.9132375 NA NA NA NA NA
## 0.25 0.6750838 0.6750838 NA NA NA NA
## 0.5 0.6294375 0.6750838 0.8742932 NA NA NA
## 1 0.6289937 0.5917343 0.8742932 0.6750838 NA NA
## 2.5 0.5917343 0.6289937 0.6750838 0.9490597 0.9490597 NA
## 5 0.2532105 0.2532105 0.6289937 0.6289937 0.6750838 0.6750838
## acet_after 0.6750838 0.6750838 0.8742932 0.8742932 0.6750838 0.9132375
## 5 acet_after
## acet_before NA NA
## 0.1 NA NA
## 0.25 NA NA
## 0.5 NA NA
## 1 NA NA
## 2.5 NA NA
## 5 NA NA
## acet_after 0.6294375 NAKW$EO4891$PvaluesMatrix## acet_before 0.1 0.25 0.5 1 2.5
## acet_before NA NA NA NA NA NA
## 0.1 0.9868483 NA NA NA NA NA
## 0.25 0.9250712 0.9868483 NA NA NA NA
## 0.5 0.6312621 0.6312621 0.6324156 NA NA NA
## 1 0.9868483 0.9868483 0.9868483 0.6312621 NA NA
## 2.5 0.6312621 0.6312621 0.6324156 0.9868483 0.6312621 NA
## 5 0.6312621 0.6891247 0.9250712 1.0000000 0.6312621 0.9868483
## acet_after 1.0000000 0.9868483 0.9868483 0.6312621 0.9868483 0.6312621
## 5 acet_after
## acet_before NA NA
## 0.1 NA NA
## 0.25 NA NA
## 0.5 NA NA
## 1 NA NA
## 2.5 NA NA
## 5 NA NA
## acet_after 0.6312621 NAKW$EO4899$PvaluesMatrix## acet_before 0.1 0.25 0.5 1 2.5
## acet_before NA NA NA NA NA NA
## 0.1 0.7298385 NA NA NA NA NA
## 0.25 0.6324156 0.6324156 NA NA NA NA
## 0.5 0.3219956 0.3219956 0.6891247 NA NA NA
## 1 0.3287615 0.3219956 0.7014248 0.7927465 NA NA
## 2.5 0.6324156 0.3219956 0.7298385 0.7298385 0.7805921 NA
## 5 0.7298385 0.7014248 0.7298385 0.6324156 0.6324156 0.6324156
## acet_after 0.6324156 0.6324156 0.7298385 0.7553550 0.7805921 0.7553550
## 5 acet_after
## acet_before NA NA
## 0.1 NA NA
## 0.25 NA NA
## 0.5 NA NA
## 1 NA NA
## 2.5 NA NA
## 5 NA NA
## acet_after 0.7298385 NAKW$EO5068$PvaluesMatrix## acet_before 0.1 0.25 0.5 1 2.5
## acet_before NA NA NA NA NA NA
## 0.1 0.8834907 NA NA NA NA NA
## 0.25 0.5393798 0.5393798 NA NA NA NA
## 0.5 0.5393798 0.8946258 0.5393798 NA NA NA
## 1 0.5393798 0.6770675 0.5393798 0.8834907 NA NA
## 2.5 0.5393798 0.5504130 0.5504130 0.6770675 0.8946258 NA
## 5 0.8946258 0.8834907 0.5393798 0.5504130 0.5504130 0.5393798
## acet_after 0.5393798 0.5393798 0.8946258 0.5393798 0.5917461 0.5393798
## 5 acet_after
## acet_before NA NA
## 0.1 NA NA
## 0.25 NA NA
## 0.5 NA NA
## 1 NA NA
## 2.5 NA NA
## 5 NA NA
## acet_after 0.5393798 NAKW$EO5663$PvaluesMatrix## acet_before 0.1 0.25 0.5 1 2.5
## acet_before NA NA NA NA NA NA
## 0.1 0.5917343 NA NA NA NA NA
## 0.25 0.5917343 0.8601193 NA NA NA NA
## 0.5 0.5917343 0.9490597 0.9132375 NA NA NA
## 1 0.7936318 0.9132375 0.9490597 0.8601193 NA NA
## 2.5 0.8601193 0.9132375 0.9132375 0.9132375 0.9132375 NA
## 5 0.9132375 0.9132375 0.8601193 0.9132375 0.9132375 0.9132375
## acet_after 0.5917343 0.9132375 0.9132375 0.9132375 0.9132375 0.8601193
## 5 acet_after
## acet_before NA NA
## 0.1 NA NA
## 0.25 NA NA
## 0.5 NA NA
## 1 NA NA
## 2.5 NA NA
## 5 NA NA
## acet_after 0.8601193 NAKW$EO5714$PvaluesMatrix## acet_before 0.1 0.25 0.5 1 2.5
## acet_before NA NA NA NA NA NA
## 0.1 0.4574500 NA NA NA NA NA
## 0.25 0.1960318 0.5256285 NA NA NA NA
## 0.5 0.2223453 0.5256285 0.7332873 NA NA NA
## 1 0.1960318 0.3429161 0.5178863 0.5256285 NA NA
## 2.5 0.3968159 0.6429369 0.7332873 0.8070399 0.8794139 NA
## 5 0.1960318 0.1960318 0.3429161 0.1960318 0.5178863 0.5178863
## acet_after 0.3429161 0.6429369 0.7332873 0.7194712 0.6429369 0.9490597
## 5 acet_after
## acet_before NA NA
## 0.1 NA NA
## 0.25 NA NA
## 0.5 NA NA
## 1 NA NA
## 2.5 NA NA
## 5 NA NA
## acet_after 0.5178863 NAKW$EO5749$PvaluesMatrix## acet_before 0.1 0.25 0.5 1 2.5
## acet_before NA NA NA NA NA NA
## 0.1 0.9441937 NA NA NA NA NA
## 0.25 0.9441937 0.9580911 NA NA NA NA
## 0.5 0.9441937 0.9441937 0.9502984 NA NA NA
## 1 0.9441937 0.9502984 0.9502984 0.9502984 NA NA
## 2.5 0.9441937 0.9441937 0.9441937 0.9502984 0.9502984 NA
## 5 0.8669914 0.8669914 0.9441937 0.9441937 0.9441937 0.9441937
## acet_after 0.8669914 0.9441937 0.9441937 0.9441937 0.9441937 0.9441937
## 5 acet_after
## acet_before NA NA
## 0.1 NA NA
## 0.25 NA NA
## 0.5 NA NA
## 1 NA NA
## 2.5 NA NA
## 5 NA NA
## acet_after 0.9502984 NAKW$EO5763$PvaluesMatrix## acet_before 0.1 0.25 0.5 1 2.5
## acet_before NA NA NA NA NA NA
## 0.1 0.9775137 NA NA NA NA NA
## 0.25 0.9775137 1.0000000 NA NA NA NA
## 0.5 0.9775137 0.9775137 0.9775137 NA NA NA
## 1 0.9775137 0.9775137 0.9775137 1.0000000 NA NA
## 2.5 0.9775137 0.9775137 0.9775137 0.9775137 1.0000000 NA
## 5 0.9775137 0.9775137 0.9775137 0.9775137 0.9775137 0.9775137
## acet_after 0.9775137 0.9775137 0.9775137 0.9775137 0.9775137 0.9775137
## 5 acet_after
## acet_before NA NA
## 0.1 NA NA
## 0.25 NA NA
## 0.5 NA NA
## 1 NA NA
## 2.5 NA NA
## 5 NA NA
## acet_after 0.9775137 NAPost hoc Wilcoxon test
Wilcoxon <- lapply(split(data, data$EO), function(i){
pairwise.wilcox.test(i$response, g = i$stimuli, p.adjust.method = "BH")
})
Wilcoxon$EO4309$p.value## acet_before 0.1 0.25 0.5 1 2.5 5
## 0.1 1 NA NA NA NA NA NA
## 0.25 1 1 NA NA NA NA NA
## 0.5 1 1 1 NA NA NA NA
## 1 1 1 1 1 NA NA NA
## 2.5 1 1 1 1 1 NA NA
## 5 1 1 1 1 1 1 NA
## acet_after 1 1 1 1 1 1 1Wilcoxon$EO4444$p.value## acet_before 0.1 0.25 0.5 1 2.5
## 0.1 0.9708625 NA NA NA NA NA
## 0.25 0.7489510 0.7489510 NA NA NA NA
## 0.5 0.7261072 0.7489510 0.9389083 NA NA NA
## 1 0.7261072 0.6812354 0.9389083 0.7489510 NA NA
## 2.5 0.6812354 0.7261072 0.7489510 1.0000000 1.000000 NA
## 5 0.2447552 0.2447552 0.7261072 0.7261072 0.748951 0.7489510
## acet_after 0.7489510 0.7489510 0.9389083 0.9389083 0.748951 0.9708625
## 5
## 0.1 NA
## 0.25 NA
## 0.5 NA
## 1 NA
## 2.5 NA
## 5 NA
## acet_after 0.7261072Wilcoxon$EO4891$p.value## acet_before 0.1 0.25 0.5 1 2.5
## 0.1 1.0000000 NA NA NA NA NA
## 0.25 1.0000000 1.0000000 NA NA NA NA
## 0.5 0.7069034 0.7069034 0.7069034 NA NA NA
## 1 1.0000000 1.0000000 1.0000000 0.7069034 NA NA
## 2.5 0.7069034 0.7069034 0.7069034 1.0000000 0.7069034 NA
## 5 0.7069034 0.7645688 1.0000000 1.0000000 0.7069034 1.0000000
## acet_after 1.0000000 1.0000000 1.0000000 0.7069034 1.0000000 0.7069034
## 5
## 0.1 NA
## 0.25 NA
## 0.5 NA
## 1 NA
## 2.5 NA
## 5 NA
## acet_after 0.7069034Wilcoxon$EO4899$p.value## acet_before 0.1 0.25 0.5 1 2.5
## 0.1 0.7856762 NA NA NA NA NA
## 0.25 0.7069034 0.7069034 NA NA NA NA
## 0.5 0.3491841 0.3491841 0.7645688 NA NA NA
## 1 0.3637607 0.3491841 0.7731546 0.8335142 NA NA
## 2.5 0.7069034 0.3491841 0.7856762 0.7856762 0.8280181 NA
## 5 0.7856762 0.7731546 0.7856762 0.7069034 0.7069034 0.7069034
## acet_after 0.7069034 0.7069034 0.7856762 0.8074095 0.8280181 0.8074095
## 5
## 0.1 NA
## 0.25 NA
## 0.5 NA
## 1 NA
## 2.5 NA
## 5 NA
## acet_after 0.7856762Wilcoxon$EO5068$p.value## acet_before 0.1 0.25 0.5 1 2.5
## 0.1 0.9288701 NA NA NA NA NA
## 0.25 0.5946524 0.5946524 NA NA NA NA
## 0.5 0.5946524 0.9314274 0.5946524 NA NA NA
## 1 0.5946524 0.7276018 0.5946524 0.9288701 NA NA
## 2.5 0.5946524 0.6012341 0.6012341 0.7276018 0.9314274 NA
## 5 0.9314274 0.9288701 0.5946524 0.6012341 0.6012341 0.5946524
## acet_after 0.5946524 0.5946524 0.9314274 0.5946524 0.6429410 0.5946524
## 5
## 0.1 NA
## 0.25 NA
## 0.5 NA
## 1 NA
## 2.5 NA
## 5 NA
## acet_after 0.5946524Wilcoxon$EO5663$p.value## acet_before 0.1 0.25 0.5 1 2.5
## 0.1 0.6812354 NA NA NA NA NA
## 0.25 0.6812354 0.9708625 NA NA NA NA
## 0.5 0.6812354 1.0000000 0.9708625 NA NA NA
## 1 0.9235431 0.9708625 1.0000000 0.9708625 NA NA
## 2.5 0.9708625 0.9708625 0.9708625 0.9708625 0.9708625 NA
## 5 0.9708625 0.9708625 0.9708625 0.9708625 0.9708625 0.9708625
## acet_after 0.6812354 0.9708625 0.9708625 0.9708625 0.9708625 0.9708625
## 5
## 0.1 NA
## 0.25 NA
## 0.5 NA
## 1 NA
## 2.5 NA
## 5 NA
## acet_after 0.9708625Wilcoxon$EO5714$p.value## acet_before 0.1 0.25 0.5 1 2.5
## 0.1 0.5310447 NA NA NA NA NA
## 0.25 0.2121212 0.5955711 NA NA NA NA
## 0.5 0.2474747 0.5955711 0.7956177 NA NA NA
## 1 0.2121212 0.3988604 0.5928516 0.5955711 NA NA
## 2.5 0.4617716 0.7132867 0.7956177 0.8666846 0.9349046 NA
## 5 0.2121212 0.2121212 0.3988604 0.2121212 0.5928516 0.5928516
## acet_after 0.3988604 0.7132867 0.7956177 0.7891502 0.7132867 1.0000000
## 5
## 0.1 NA
## 0.25 NA
## 0.5 NA
## 1 NA
## 2.5 NA
## 5 NA
## acet_after 0.5928516Wilcoxon$EO5749$p.value## acet_before 0.1 0.25 0.5 1 2.5
## 0.1 0.9946531 NA NA NA NA NA
## 0.25 0.9946531 1.0000000 NA NA NA NA
## 0.5 0.9946531 0.9946531 0.9946531 NA NA NA
## 1 0.9946531 0.9946531 0.9946531 0.9946531 NA NA
## 2.5 0.9946531 0.9946531 0.9946531 0.9946531 0.9946531 NA
## 5 0.9790210 0.9790210 0.9946531 0.9946531 0.9946531 0.9946531
## acet_after 0.9790210 0.9946531 0.9946531 0.9946531 0.9946531 0.9946531
## 5
## 0.1 NA
## 0.25 NA
## 0.5 NA
## 1 NA
## 2.5 NA
## 5 NA
## acet_after 0.9946531Wilcoxon$EO5763$p.value## acet_before 0.1 0.25 0.5 1 2.5 5
## 0.1 1 NA NA NA NA NA NA
## 0.25 1 1 NA NA NA NA NA
## 0.5 1 1 1 NA NA NA NA
## 1 1 1 1 1 NA NA NA
## 2.5 1 1 1 1 1 NA NA
## 5 1 1 1 1 1 1 NA
## acet_after 1 1 1 1 1 1 1EAG plot
for some reason, i cannot add the leg ID to the hovering text in the box plot, but i can do it in the dot-plot. so in order to detect outliers, we can detect them in the first plot (the box plot), then identify their ID leg in the dot plot:)
dot <- ggplot(data, aes(x = stimuli, y = response, colour = EO, leg=leg)) +
geom_point() +
facet_wrap( ~ EO) +
theme_linedraw() +
ggtitle("Varroa foreleg response to different essential oils") +
xlab("Stimuli amount (microgram)") +
ylab("Normalized response (%)") +
theme(axis.text.x=element_text(angle=45, hjust=1)) +
stat_summary(aes(group=EO), fun=mean, geom="line", colour="black")
ggplotly(dot, tooltip = c("leg","response"))